Critical behavior and synchronization of discrete stochastic phase coupled oscillators 
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Synchronization of stocliastic phase-coupled osciUators is known to occur but difficult to character- 
ize because sufficiently complete analytic work is not yet within our reach, and thorough numerical 
description usually defies all resources. We present a discrete model that is sufficiently simple to 
be characterized in meaningful detail. In the mean field limit, the model exhibits a supercritical 
Hopf bifurcation and global oscillatory behavior as coupling crosses a critical value. When coupling 
between units is strictly local, the model undergoes a continuous phase transition which we charac- 
terize numerically using finite-size scaling analysis. In particular, we explicitly rule out multistability 
and show that that the onset of global synchrony is marked by signatures of the XY universality 
class. Our numerical results cover dimensions d = 2, 3, 4, and 5 and lead to the appropriate XY 
classical exponents /3 and v, a lower critical dimension die ~ 2, and an upper critical dimension 
due = 4. 

PACS numbers: 64.60.Ht, 05.45.Xt, 89.75.-k 



INTRODUCTION 



The role of dissipative structures and self-organization 
in systems far from equilibrium in the description of real 
and observable physical phenomena has been undisputed 
since the experiments with the Belusov-Zhabotinsky re- 
actions in the early 1960's. The breaking of time transla- 
tional symmetry has since become a central and typical 
theme in the analysis of nonlinear noncquilibrium sys- 
tems. It is somewhat surprising that in the later stud- 
ies of spatially distributed systems, most of the interest 
shifted to pattern forming instabilities, and little atten- 
tion was devoted to the phenomenon of bulk oscillation 
and the required spatial frequency and phase synchro- 
nization, especially in view of the intense interest gen- 
erated in the scientific and even broader community by 
the emergence of phase synchronization in populations of 
globally coupled phase oscillators 1]. The synchronous 
firing of firefiies is one of the most visible and spectac- 
ular examples of phase synchronization. Because intrin- 
sically oscillating units with slightly different eigenfre- 
quencies underlie the macroscopic behavior of an exten- 
sive range of biological, chemical, and physical systems, 
a great deal of literature has focused on the mathemat- 
ical principles governing the competition between indi- 
vidual oscillatory tendencies and synchronous coopera- 
tion H, H, Q . While most studies have focused on glob- 
ally coupled units, leading to a mature understanding of 
the mean field behavior of several models, relatively lit- 
tle work has examined populations of oscillators in the 
locally coupled regime 0- In fact, models of lo- 

cally coupled oscillators typically involve a prohibitively 
large collection of interdependent nonlinear differential 
equations, thus preventing any extensive characteriza- 



tion of the phase transition to phase synchrony. Fur- 
ther inclusion of stochastic fluctuations in such models 
typically renders them computationally and analytically 
intractable for even a modest number of units. As a re- 
sult, the description of emergent synchrony has largely 
been limited to small-scale and/or globally-coupled de- 
terministic systems Q , despite the fact that the dy- 
namics of the physical systems in question likely reflect a 
combination of finite-range forces and stochasticity. Two 
recent studies by Risler et al. 0, Inj represent notable 
exceptions to this trend. Using an elegant renormaliza- 
tion group approach, they provide analytical evidence 
that identical locally-coupled noisy oscillators belong to 
the XY universality class, though to date there had been 
no empirical verification, numerical or otherwise, of their 
predictions. 

The difficulty with existing models of locally coupled 
oscillators is that each is typically described by a non- 
linear differential equation, and the resulting systems of 
coupled equations are computationally extravagant, es- 
pecially when stochastic components are also included. 
Here we introduce a far more tractable model consisting 




FIG. I: Single three state unit with generic transition rates 
9- 
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of identical and discrete phase-coupled oscillators whose 
simple structure renders it amenable to extensive numer- 
ical study. The use of such minimal models is common 
in statistical physics, where microscopic details can often 
be disregarded in favor of phenomenological macroscopic 
variables. As Landau theory reminds us, macroscop- 
ically observable changes (those that occur on length and 
time scales encompassing a magnificently large number 
of degrees of freedom) occur without reference to mi- 
croscopic specifics. In a sense, the distinguishing fea- 
tures of even highly diverse systems become irrelevent 
for the description of cooperative behavior at the level 
of a phase transition; instead, the underlying statisti- 
cal similarities give rise to classes of universal behavior 
whose members differ greatly at the microscopic level. 
In the spirit of this universality, simple toy models have 
been devised in hopes of capturing the essential quali- 
tative features of phase transitions without concern for 
the microscopic structure of the problem. With this 
in mind, we construct the simplest model that exhibits 
global phase synchrony and contains the physical ingre- 
dients listed above, namely, stochastic variation within 
individual units and short-ranged interactions [l^ . The 
simplicity of the model allows for relatively fast numer- 
ical simulation and thus an extensive description of the 
phase transition in question. An abbreviated version of 
our principal results can be found in 13]. There we char- 
acterized the universality class of the transition, includ- 
ing the critical exponents and the lower and upper critical 
dimensions. Here we present considerably more detail as 
well as additional results to support our characterization. 

The organization of the paper is as follows. In Sec. UTI 
we introduce our description of a single unit as well as 
the coupling scheme between units. Section IlIII presents 
the linear stability analysis of the mean field limit, and 
Sec. IIVI contains the finite-size scaling analysis that un- 
veils the critical behavior of the locally-coupled model. 
We conclude with a summary in SeclVl 



II. THREE STATE MODEL 

Our starting point is a three-state unit governed 
by transition rates g, as shown in Fig. ^ We interpret 
the state designation as a generalized (discrete) phase, 
and the transitions between states, which we construct 
to be unidirectional, as a phase change and thus an oscil- 
lation of sorts. The probability of going from the current 
state i to state i -I- 1 in an infinitesimal time dt is gdt, 
with i = 1,2,3 modulo 3. For a single unit, g is simply a 
constant that sets the oscillator's intrinsic frequency; for 
many units coupled together, we will allow g to depend 
on the neighboring units in the spatial grid, thereby cou- 
pling neighboring phases. The choice of coupling, speci- 
fied below, is not unique. 

For a single unit we write the linear evolution equation 



d_ 
di 



P{t) = MP{t) 



where 



(1) 



(2) 



Piit) is the probability of being in state i at time i, and 



M 




(3) 



The system clearly reaches a steady state for P* — Pj* = 
PI = 1/3. The transitions 1^2, 2^3, 3^1 occur 
with a rough periodicity determined by g. The time evo- 
lution of our simple model thus qualitatively resembles 
that of the discretized phase of a generic noisy oscillator. 

We are interested in the behavior that emerges when 
individual units are coupled to one another by allowing 
the transition probability of a given unit to depend on 
the states of the unit's nearest neighbors in the spatial 
grid. The phase at a given site is compared with those of 
its neighbors, and the phase of the given site is adjusted 
so as to facilitate phase coherence. The expectation is 
to capture the physical nature of synchronization. It is 
further expected that within certain restrictions (e.g., the 
coupling must surely be nonlinear), the specific nature of 
the coupling is not important (in other words, we expect 
universality) so long as we ultimately observe a transition 
to global synchrony at some finite value of the coupling 
parameter. We settle upon a particular exponential form 
below. As we shall see, linear stability analysis for this 
choice confirms the existence of a Hopf bifurcation in the 
mean field limit. 

More concretely, we specify that each unit may tran- 
sition to the state ahead or remain in its current state 
depending on the states of its nearest neighbors. For 
unit /i, which we take to be in state «, we choose the 
transition rate to state j as follows: 



5ij = gexp 



2d 



i+l, 



(4) 



where a is the coupling parameter, 6 is the Kronecker 
delta, TVfe is the number of nearest neighbors in state fc, 
and 2c? is the total number of nearest neighbors in d di- 
mensions. The transition rate is thus determined by the 
number of nearest neighbors of unit /i that are one state 
ahead and in the same state as unit /i. Table I shows the 
explicit transition rates in one dimension. While these 
rates are somewhat distorted by their assumed indepen- 
dence of the number of nearest neighbors in state i — 1 
(e.g. in one dimension the transition rate from state i to 
state z -I- 1 is the same if both nearest neighbors are in 
state i — 1 and if one is in state i and the other in i -I- 1), 
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the form Q is simplified by this assmumption and, as we 
shall see, does lead to synchronization. Note also that an 
equally simple model might posit a coupling which de- 
pends on Ni^i, the number of units 'behind' the unit 
in question, rather than A'^^+i, or a more complex model 
could be constructed that depends on both. We settle 
on our choice Q because the phase transition we seek 
occurs for a smaller value of the coupling constant a, and 
therefore numerical simulations can be run with larger 
time steps. We stress again that universality suggests 
that such microscopic details should not substantially al- 
ter the qualitative picture of the phase transition as long 
as the coupling is sufficiently nonlinear and favors syn- 
crhonization. 
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III. MEAN-FIELD THEORY 

To test for the emergence of global synchrony, we first 
consider a mean field version of the model, that is, one 
where each unit is coupled to all other units. In the large 
A'^ limit with all-to-all coupling we write 



gij = g exp [a{Pj - Pi)] 6. 



(5) 



where Pk is the (ensemble) probability of being in state 
k. Note that with all-to-all coupling gij does not depend 
on the location of the unit within the lattice. Note also 
that there is an inherent assumption that we can replace 
Nk/N with Pk, that is, we are assuming that N, the total 
number of units, is large enough that Nk/N serves as a 
good estimation of the ensemble probability Pk- With 
this simplification we arrive at an equation for the mean 
field P: 



^^P{t)^M[P{t)]P{t), 



where 



M[P{t)] = 



' -g\2 .931 

.912 -323 
\ 323 -531, 



(6) 



(7) 
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FIG. 2: Simulations with 5000 globally coupled units (bottom 
panel) agree well with the numerical solution of the mean field 
equations (top panel). As predicted by linearization, a Hopf 
bifurcation occurs near a = 1.5. 



We have explicitly noted the dependence of M on P{t) 
since each of the matrix elements gij depends on the 
evolving probabilities. Equation (|SJ) is thus a highly non- 
linear equation. 

The normalization condition Pi + -P2 + ^3 — 1 allows us 
to eliminate P3 and obtain a closed set of equations for Pi 
and P2. Wc can further characterize the mean field solu- 
tions using standard linear stability analysis. Specifically, 
we linearize about the fixed point {Pi,P2) = (V^i 1/3) 
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TABLE I: Transition rates in one dimension. 



FIG. 3: Absence of synchronization in 2D. Top left: a — 1.5. 
Top right: a = 2.5. Bottom left: a — 3.5. Bottom right: 
a = 4.5. L = 100 for all plots. Even for very large values 
of the coupling, highly synchronous oscillatory behavior is 
not present. As discussed in the text and shown in the next 
figure, the intermittent oscillations apparent for high values 
of a result from finite size effects. 
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and obtain the Jacobian J evaluated at {P* , P2)' 

V .9 - .9/ 

The eigenvalues of J characterize the fixed point, 
they are given by: 

A± = |(2a - 3±iV3). 



(8) 



and 



(9) 



Both cross the imaginary axis at a = 1.5, indicative of 
a Hopf birfurcation at this value. Hence, as a increases 
the mean field undergoes a qualitative change at a = 1.5 
from disorder (Pi = P2 = P3) to global oscillations, and 
the desired global synchrony emerges. 

The predictions of the linearization can be verified by 
numerically solving the mean field equations ©. In 
turn, these solutions agree well with direct simulations 
of the multiple unit model characterized by Eq. if we 
consider all-to-all coupling rather than merely nearest- 
neighbor coupling (Fig. [21 . As such, the mean field 
equations accurately capture the behavior of the nearest 
neighbor model in the high (spatial) dimensional limit. 

Furthermore, the Hopf bifurcation seen in our mean 
field model can be shown to be supercritical. Such an 
analytical argument is formally related to the structure 
of the normal form for the Hopf bifurcation. Practically 
speaking, one must consider the sign of the first Lya- 
punov coefficient at the bifurcation point (oc = 1.5). 
Following (l5l |. we transform our two dimensional non- 
linear equation © to a single equation for the complex 
variable z valid for small a — a ^ Qc- The form of the 
equation is given by 



z — \{a)z + f{z, z^, a), 



(10) 
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FIG. 4: Log- log plot of r vs for d — 2. The order parame- 
ter r tends to as system size increases, verifying the absence 
of a transition in two dimensions. Even for large values of 
the coupling, synchronous oscillations die away in the limit of 
infinite system size. 



where f{z,z''^,a) is an 0(|zp) smooth function of z, z\ 
and a, and Xia) is an a-dependent eigenvalue of the lin- 
earized Jacobian (|SJ) given above. We achieve such a 
transformation by first finding complex eigenvectors p 
and q given by 



J(0)g = A(0)g,J(0)^p = A(0)tp, 



(11) 



with J(0) the Jacobian evaluated at a = Oc = 1.5. One 
then normalizes (p, q) , where brackets in this context rep- 
resent the standard complex scalar product. An equation 
for z of the desired form 11U|) is formally attained at a = 
as 

z-A(0)z+(p,F(z9 + ;^t^^0)), (12) 

where F(x^ a) is related to our original dynamical system, 
i.e., 

^P(t) = J{a)P{t)^F\P{t\a\- (13) 

From this, we may obtain the first Lyapunov coefficeint 
L\ as 

= ;T-2-Re(«/2o/ii + Wo/21), (14) 
with jij given by the formal Taylor expansion of /, 
/(z,zt,o) = (p,F(z<z + zV^O))= ^ _L/,,^fe(^ty. 

/c+/>2 

(15) 

An explicit calculation for our mean field dynamical sys- 
tem reveals that Li < 0, indicative of a supercritical Hopf 
bifurcation to a unique, stable limit cycle as a eclipses Oc. 

In what follows, we characterize the breakdown of 
the mean field description as spatial dimension is de- 
creased, and characterize the phase transitions observed 
with nearest-neighbor coupling. 



IV. CRITICAL BEHAVIOR OF THE LOCALLY 
COUPLED MODEL 

With a firm understanding of the mean field model, we 
now follow with a study of the locally coupled model. We 
perform simulations in continuous time on d-dimensional 
cubic lattices of various sizes. For all simulations, we im- 
plement periodic boundary conditions. Time steps dt are 
taken to be 10 to 100 times smaller than the fastest pos- 
sible local average transition rate, that is, dt ^ e"" (we 
set g — \ \vl our simulations). This estimate is actually 
quite conservative, particularly because the fastest possi- 
ble transition corresponds to a single unit in state i with 
all Id nearest neighbors in state z -I- 1, a scenario which 
certainly does not dominate the macroscopic dynamics. 
We have ascertained that differences between these sim- 
ulations and others run at much smaller time steps (500 
to 1000 times smaller than e~°) are very small. All 
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FIG. 5: Snapshots of the system in d = 2 are shown for a = 1, 
a — 2.5, a = 4.5, and a = 6.5. Upon close inspection, one 
can discern vortex-like structures, particularly for the higher 
values of a. The three colors represent units in the three 
possible states. 



simulations were run until an apparent steady state was 
reached. Furthermore, we start all simulations from ran- 
dom initial conditions, and we calculate statistics based 
on 100 independent trials. Although the simplicity of 
the model allows for efficient numerical simulation, our 
results nevertheless represent a modest computational 
achievement; simulations required approximately 5 weeks 
on a 28-node dual processor cluster. 

To characterize the emergence of phase synchrony, we 
introduce the order parameter Q 



R 



1 ^ 

-IF 



I 



(16) 



Here is a discrete phase, taken to be 27r(A; — l)/3 for 
state k G {1,2,3} at site j. The brackets represent an 
average over time in the steady state and an average over 
all independent trials. A nonzero value of r in the ther- 
modynamic limit signifies the presence of synchrony. We 
also make use of the corresponding generalized suscepti- 
bility 



(17) 



We begin by considering the model in two spatial di- 
mensions. Here, as shown in Fig. |21 we do not see the 
emergence of global oscillatory behavior. Instead, we ob- 
serve intermittent oscillations (for very large values of a) 
that decrease drastically with increasing system size. In 
fact, as depicted in Fig.^ r approaches zero in the ther- 
modynamic limit, even for very large values of a. We 
conclude that the phase transition to synchrony cannot 
occur for d = 2. Interestingly, snapshots of the sytem re- 
veal increased spatial clustering as a is increased as well 



as the presence of defect structures, perhaps indicative of 
Kosterlitz-Thouless-type phenomena (Fig. O Fur- 
ther studies along these lines are underway. 

In contrast to the d = 2 case, which serves as the lower 
critical dimension, a clear thermodynamic-like phase 
transition occurs in three spatial dimensions. We see the 
emergence of global oscillatory behavior, which persists 
in the limit of large system size, as a increases past a 
critical value Uc (Figs. El and [Tj). This is consistent with 
the predictions of the mean field theory. For a < ac, r 
approaches zero as system size is increased, and a dis- 
ordered phase persists in the thermodynamic limit. As 
expected, for a > Uc the steady state dynamics of Pi and 
P2 exhibit smooth temporal oscillations (see the lower 
insets in Fig. 0) similar to the mean field case beyond 
the Hopf bifurcation point. In addition, Fig. [7| shows 
the behavior of the order paramater as a is increased for 
the largest system studied {L — 80); the upper left inset 
shows the peak in x at a = 2.345 ± 0.005, thus provid- 
ing an estimate of the critical point ac where fluctuations 
are largest. Strictly speaking, we must extrapolate this 
peak to obtain a result in the limit of infinite system size, 
but we see no change as system size is increased beyond 
L = 40, indicating that finite size effects are small in sys- 
tems beyond this size. At any rate, such finite size effects 
are within the range of our estimation. We tried to apply 
the Binder cumulant crossing method 16] for determin- 
ing flc more precisely, but residual finite size effects and 
statistical uncertainties in the data prevent us from de- 
termining the crossing point with more precision than 
that stated above. In any case, we are only interested 
in determining the critical point with sufficient accuracy 
to determine the universality class of the transition. For 
this, as we show below, our current estimation suffices 
in three dimensions as well as in higher dimensions. In 
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FIG. 6: Log-log plots of r vs for d = 3. For a > Oc, 
the order parameter r approaches a finite value, even as the 
system size increases indefinitely. For a < Uc, r approaches 
zero in the thermodynamic limit. 
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FIG. 7: Onset of synchronization in d = 3. Global oscillatory 
behavior emerges as a is increased beyond flc, as indicated by 
the increasing value of r. The system size is L = 80. Upper 
left inset: Fluctuations peak near the critical point, giving 
an estimation of ac = 2.345 ± 0.005. Right insets: Pi and 
P2 undergo smooth temporal oscillations for large a (upper 
right), while a lower value of a decreases temporal coherence 
(lower left). 



addition, we note that the transiition to synchrony ap- 
pears to be a smooth, second order phase transition. To 
rule out potential multistability (and thus a discontinu- 
ous first order transition), we show histograms of r for 
d = 4 given over all independent trials in Fig. |H1 The his- 
tograms show no evidence whatsoever of multiple peaks 
beyond the statistical fluctuations expected for the rela- 
tively small sample size, and thus we can safely rule out a 
discontinuous transition, in agreement with the findings 
of the mean field analysis. Similarly peaked histograms 
are found in d = 3 (less sharply peaked but distinctly 
unimodal) and d — 5 (more sharply peaked). 

To further characterize this transition, we use a sys- 
tematic finite size scaling analysis. We start by assuming 
the standard form for r in a finite system, 



r^L--F[{a~ac)L^], 



(18) 



where F{x) is a scaling function that approaches a con- 
stant as a; — > 0. This ansatz suggests that near the crit- 
ical point we should plot rL^ vs. {-^ — and data 
from different system sizes should collapse onto a sin- 
gle curve. To test our numerical data against different 
universality classes we choose the appropriate critical ex- 
ponents for each, recognizing that there are variations 
in the reported values of these exponents. For the XY 
universality class we use the exponents reported in 
(/3 = 0.34 and v = 0.66). For the Ising exponents we use 
those given in [l8| (/? = 0.31 and i' = 0.64). In Fig.El we 
see quite convincingly a collapse when exponents from 
the XY class are used. For comparison, we also show the 
data collapse when 3D Ising exponents are used (Fig. llO|l . 
Our data suggests that the model falls within the XY 
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FIG. 8: Lack of multistability in d = 4: Histograms over all 
independent trials show only single peaks of varying widths, 
consistent with the expectations for a second order phase 
transition. 
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FIG. 9: Exponents in d = 3: Data collapse of rL " vs {a/a^ — 
1)1/'' . With ttc = 2.345, we show the data collapse using the 
theoretical XY exponents in 3D. The collapse is excellent, 
suggesting that the model is in the XY universality class. 
The insets show a closer view near the critical point. 



Universality class, though the very small differences be- 
tween XY and Ising exponents makes it impossible to 
entirely rule out Ising-like behavior. We should point 
out that while some reported values of the Ising criti- 
cal exponents differ from the XY values by more than 
those used above, others differ by less (see [l^ for an ex- 
haustive collection of estimates). Note that this scaling 
procedure was attempted for many values a^. within the 
stated range of accuracy. In all cases where a distinc- 
tion could be made, the XY exponents provided a better 
collapse than the corresponding Ising exponents. 

To complete the analogy with the equilibrium phase 
transition, we explore spatial correlations in c? = 3. 
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FIG. 10: Exponents in d = 3. Data collapse of rL " vs (a/flc — 
1)L^. With ac = 2.345, we show the data collapse using 
theoretical Ising 3D exponents. The collapse is reasonable 
good, but still poor compared with that seen with exponents 
from the XY class. Insets show a closer view near the critical 
point. 



FIG. 11: Spatial correlations in d = 3. As a approaches the 
critical value Oc, evidence of long range correlations devel- 
ops, indicative of a diverging correlation length at the critical 
point. The lower inset shows the power law decay of the cor- 
relation function at the critical point, while the upper inset 
shows that the correlation function decays exponentially far 
from the critical point. 



Specifically, we calculate C(Z), the spatial correlation 
function, given by 

N 3 

^(0 = (E E ^'^i) ^-^^i+O) - r\ (19) 

j=l n=l 

Here again indicates the discrete phase of the oscilla- 
tor at site j, and Z„ denotes the Cartesian components in 
the X, y, and z directions at distance I from site j. The 
correlation function depends only on this distance. As 
seen in Fig. 1111 correlations develop for values of a near 
the critical point, while this correlation is absent away 
from Oc. The functional form of C{1) as a approaches 
flc is similar to that seen in equilibrium transitions. In- 
deed, the lower inset is at the critical point (a — 2.345) 
and explicitly shows power law decay of the correlation 
function. The upper inset is far from the critical point 
(a = 1.8) and shows exponential decay. 

In four spatial dimensions we also see a transition to 
synchrony characterized by large fluctuations at the crit- 
ical point. Here we estimate the transition coupling to 
be Oc = 1.900 ± 0.025 by again considering the peak 
in X (see Figs. WI\ and [T^ . Because we expect d = 4 
to be the upper critical dimension in accordance with 
XY/Ising behavior, we anticipate a slight breakdown of 
the scaling relation (|18|l . An alternate scaling ansatz 
valid at due is given by (|18|l with the transformation 
L — > ln{L)L^/^ U^- ^ priori it is not clear how strongly 
(|18|l will be violated in d = 4, nor is it clear that the 
modified ansatz will better serve our purposes; therefore, 
we will use both forms of scaling in testing for the mean 
field exponents v = 1/2 and (3 — 1/2. 

As shown in Fig. ^1 the data collapse is very good 
with the mean field exponents regardless of which scaling 
ansatz is used. As such, our simulations suggest that 
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a 

FIG. 12: Transition in d = 4: The behavior of the order pa- 
rameter near the transition point is shown for various system 
sizes. The inset shows the generalized susceptablitiy, X) which 
peaks at a = 1.900 ± 0.025, giving an estimation of Uc- 



d = 4 serves as the upper critical dimension; additionally, 
it appears that corrections to finite-size scaling at d = 4 
are not substantial, though a much more precise study 
would be needed to investigate such corrections in greater 
detail. 

To further support the claim that due = 4, we consider 
the case d = 5. We see a transition to synchrony which 
occurs at ae = 1.750 ± 0.015 (see Figs. [T51 and TO. As 
expected, this value for ae is considerably closer than the 
critical coupling in four dimensions to the value Oc = 
1.5 calculated by linear stability analysis in mean field 
theory. 

Finally, it is interesting to test the suggestion of Jones 
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FIG. 13: Log-log plots of r vs for d = 4. The order 
parameter r clearly approaches a finite, nonzero value for a > 
flc and approaches for a < ac. 
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peaks at a = 1.750 ± 0.015, giving an estimation of Uc- 
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FIG. 14: Exponents in d = 4: Data collapse of original 
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FIG. 16: Log-log plots of r vs in d — 5. The order 
parameter r clearly approaches a finite, nonzero value for a > 
ac and approaches for a < Uc- The value of Oc appears to 
fall between a = 1.8 and a = 1.7. 



and Young 0| that above the critical dimension, d ^ 
due, it is appropriate to modify the finite size scaling 
ansatz (|I8f) by the transformation L L'^/'^. We test 
this suggestion for d = 5. As indicated in Fig. El the 
data collapse is excellent for both the original scaling 
and the modified form of the ansatz. The collapse of 
the data with mean field exponents seems slightly better 
using the modified ansatz, though a much more precise 
study would be required to accurately capture the form 
of the modified scaling in d > d^c- In any case, our data 
suggest that the model exhibits mean field behavior in 
d = 5, verifying that d — A serves as the upper critical 
dimension. 



V. SUMMARY 

We have introduced a simple discrete model for study- 
ing phase coherence in spatially distributed populations 
of noisy coupled oscillators. This model lends itself to 
numerical study even in the case of nearest neighbor cou- 
pling because each oscillator is a simple three-state sys- 
tem rather than one of the usual continuum choices. The 
coupled system is therefore much simpler than the usual 
set of coupled nonlinear differential equations. 

A mean field treatment combined with linear stability 
analysis shows that the globally coupled model undergoes 
a Hopf bifurcation to macroscopic synchrony as the cou- 
pling parameter a is increased. We are able to determine 
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collapse of the data is quite convincing when the exact mean 
field exponents are used. 



the mean field critical coupling constant analytically. For 
locally coupled units, numerical solution of the system 
shows the emergence of a thermodynamic synchronous 
phase for d > 2, indicating that the lower critical di- 
mension is die — 2. As d is increased, the numerically 
established critical value Oc approaches that predicted 
by the mean field treatment of the model. For d = 3, we 
give strong numerical evidence that the model falls into 
the 3D XY universality class, while for d = 4 the criti- 
cal exponents are those predicted by mean field theory. 
The exponents in c? = 5 also take on the mean field val- 
ues, thus verifying that d = 4 corresponds to the upper 
critical dimension due- 



In conclusion, while nonequilibrium phase transitions 
have a much wider diversity in universality classes than 
equilibrium ones , it is remarkable that the prototype 
of a nonequilibrium transition, namely, a phase transi- 
tion that breaks the symmetry of translation in time, 
is described, at least for the critical exponents investi- 
gated in this paper, by an equilibrium universality class. 
In particular, the Mermin- Wagner theorem, stating that 
continuous symmetries can not be broken in dimension 
two or lower, appears to apply. Furthermore, the XY 
model is known to display a Kosterlitz-Thouless tran- 
sition, in which, beyond a critical temperature, vortex 
pairs can unbind into individual units creating long range 
correlations. Preliminary results indicate that a similar 
transition occurs in our model. Finally, a note of cau- 
tion concerning the discreteness of the phase is in or- 
der. We first note that microscopic models often feature 
discrete degrees of freedom. For example, our model is 
reminiscent of the triangular reaction model introduced 
by Onsager ^221 , on the basis of which he illustrated the 
concept of detailed balance as a characterization of equi- 
librium. Continuous phase models appear in a suitable 
thermodynamic limit. We stress that the breaking of 
time translational symmetry can occur independently of 
whether the phase is a discrete or continuous variable. It 
is, however, not evident whether continuous and discrete 
phase models belong to the same universality class. The 
results found here seem to support the latter thesis, but 
a renormalization calculation confirming this hypothesis 
would be welcome. 
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